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Merging supermassive black hole-black hole (BHBH) binaries produced in galaxy mergers are 
promising sources of detectable gravitational waves. If such a merger takes place in a gaseous en- 
vironment, there is a possibility of a simultaneous detection of electromagnetic and gravitational 
radiation, as the stirring, shock heating and accretion of the gas may produce variability and en- 
hancements in the electromagnetic flux. Such a simultaneous detection can provide a wealth of 
opportunities to study gravitational physics, accretion physics, and cosmology. We investigate this 
scenario by performing fully general relativistic, hydrodynamic simulations of merging, equal-mass, 
nonspinning BHBH binaries embedded in gas clouds. We evolve the metric using the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) formulation with standard moving puncture gauge conditions 
and handle the hydrodynamics via a high-resolution shock-capturing (HRSC) scheme. We consider 
both "binary Bondi accretion" in which the binary is at rest relative to the ambient gas cloud, 
as well as "binary Bondi-Hoyle-Lyttleton accretion" in which the binary moves relative to the gas 
cloud. The gas cloud is assumed to be homogeneous far from the binary and governed by a T-law 
equation of state. We vary Y between 4/3 and 5/3. For each simulation, we compute the gas flow and 
accretion rate and estimate the electromagnetic luminosity due to bremsstrahlung and synchrotron 
emission. We find evidence for significant enhancements in both the accretion rate and luminosity 
over values for a single black hole of the same mass as the binary. We estimate that this luminos- 
ity enhancement should be detectable by LSST for a 1O 6 M0 binary in a hot gas cloud of density 
n ~ 10 cm -3 and temperature T ~ 10 6 K at z = 1, reaching a maximum of L ~ 3 x 10 43 erg s _1 , 
with the emission peaking in the visible band. 

PACS numbers: 04.25.D-, 04.25. dg, 47.75. +f 



I. INTRODUCTION 

All bulge galaxies (including the Milky Way) are 
believed to contain a central supermassive black hole 
(SMBH) with a mass M between 10 4 M Q and 1O 9 M 
ELS Hi- It is also believed that galaxy mergers commonly 
lead to the formation of a massive black hole-black hole 
(BHBH) binary system in the merged remnants [3, Q . In 
the standard picture, the BHBH binary separation de- 
creases first through dynamical friction due to distant 
stellar encounters, then through gravitational slingshot 
interactions in which nearby stars are ejected at velocities 
comparable to the binary's orbital velocity, and finally 
through the emission of gravitational radiation, leading 
to coalescence of the binary Q. These low- frequency 
gravitational waves will be detectable by LISA and will 
contain a wealth of information about the inspiral. How- 
ever, it has been argued that the gaseous accretion flow 
which forms around the binary can be a source of elec- 
tromagnetic radiation as well Q . This raises the exciting 
possibility of a simultaneous detection of electromagnetic 
and gravitational waves from BHBH mergers. 

This picture is supported by a number of observed 
AGNs that may be harboring BHBH binaries. VLBI ob- 
servations of the elliptical galaxy 0402+379 have discov- 
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ered two radio sources at a projected separation of only 
7 pc. The existence of jets, as well as variability asso- 
ciated with BH activity, indicate that the system may 
be a BHBH binary [H, 0, [1(3]. Another candidate is OJ 
287, a BL Lac object whose light curve shows variabil- 
ity with a period of ~ 12 yr. It is believed that this 
may be a massive BHBH binary in which the smaller 
BH orbits with a period of 12 yr, penetrating the disk 
of t he p rimary and giving rise to the observed variabil- 
ity [HI [H [T|. The quasar SDSS 092712.65+294344 
is not believed to be a binary system, but rather a re- 
coiling BH which is the product of a binary merger. 
This suggestion is supported by a systematic shift of 
2650 km s _1 in its emission lines [III [la]- Another candi- 
date is quasar SDSS J153636. 22+044127.0, in which two 
broad line emission systems are observed, separated in 
velocity by 3500 km s _1 . This observation has been in- 
terpreted as a BHBH binary system in which each object 
has its own emission system [l6[. 

Information from a simultaneous detection of elec- 
tromagnetic and gravitational waves may be useful for 
studying fundamental aspects of gravitational physics. 
For example, in some modified gravity scenarios, the 
propagation velocity for gravitons may differ from that 
of photons [TtI [HI] • Additionally, the measurement of the 
luminosity distance from the gravitational wave signal at 
an accuracy of 1—10%, coupled with the redshift informa- 
tion from the electromagnetic detection, could serve as a 
cosmological "standard siren" of unprecedented accuracy 
(better than ~ 1%) [l9(. Such detections may also com- 
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bine accurate measurements of BH spins and masses ob- 
tained from gravitational wave signals with electromag- 
netic observations to probe BH accretion physics in great 
detail [H. 

Several mechanisms for electromagnetic emission from 
accretion disks around merging BHBH binaries have been 
proposed. In one scenario, the inner edge of the accre- 
tion disk is identified as the radius at which the viscous 
torque on the gas balances the gravitational torque from 
binary. This radius is between 1.5 and 2 times the orbital 
separation [2l|, |22j, |23, |2J] and encompasses a hollowed 
region in the disk. Late in the inspiral the BHBH binary 
decouples from the disk and coalesces. For a binary of 
mass M « 10 6 M o accreting at 10% of the Eddington 
rate, the subsequent evolution of this disk, assumed op- 
tically thick, gives rise to a source which initially peaks 
in the UV band and hardens to EUV and soft X-ray 
emission at late times 0, Additionally, there is a 

sudden change in the mass of the binary during the final 
stage of the merger, as gravitational waves carry away 
a few percent of the mass. The abrupt change in the 
potential creates waves in the disk which can grow into 
shocks and increase the luminosity of the disk in a char- 
acteristic way [H, [27], HH, giving rise to a potentially 
detectable prompt X-ray signal. Another possibility is 
that the merged BH remnant experiences a recoil veloc- 
ity which ma y, in principle, be as high as several thou- 
sand km s _1 [20| , although it is likely to be much lower 
(< 200 km/s) in most galaxy mergers [30]. This recoiling 
BH may "shake" or penetrate the disk, creating shocks 
which could give rise to a transient signal. 

Various methods have been used to model plausible 
sources of electromagnetic emission from BH mergers. In 
one approach, the dynamics of the inspiral is ignored, fo- 
cusing on the effect of BH kicks and/ or BH mass loss on 
the hydrodynamical flow [H, H3, [H, HH II H, H, 11 . 
In another approach, the behavior of the gas is modeled 
by following the motion of collisionless "particle tracers" 
on geodesies [36|. Only recently has a fully relativistic, 
dynamical simulation of a BHBH binary merger in a hy- 
drodynamic setting been performed [371 ]. 

In this paper we study BHBH binary mergers in the 
presence of ambient gas. Modeling such systems re- 
quires fully general relativistic dynamical simulations, in- 
cluding relativistic hydrodynamics. The development of 
stable algorithms to integrate Einstein's field equations 
of general relativity numerically in 3 + 1 dimensions, 
such as the BSSN formalism [33, and the general- 
ized harmonic approach [40j |. together with the identifi- 
cation of suitable gauge conditions, has enabled several 
pioneering simulations that demonstrate how to track 
the late ins pira l, pl ung e and merger of a BHBH binary 
in vacuum [4ll . |42j, |43j . More refined follow-up simula- 
tions of these strong-field, late phases, joined onto an- 
alytic, post-Newtonian calculations of the early inspiral 
epoch [441 ], are now capable of producing accurate gravi- 
tational waveforms for merging BHBH binaries with com- 
panions spanning a range of mass ratios and spins. 



With the problem of gravitational wave emission from 
a vacuum BHBH binary inspiral well in hand, it is now 
important to turn to the problem of electromagnetic 
emission from BHBH binary coalescence in an astrophys- 
ically realistic gaseous environment. 

When gas accretes onto a SMBH in the center of a 
galaxy the specific angular momentum of the gas L may, 
in many cases, be much larger than that of a circular 
orbit near the horizon, L <~ 2Mc [45[ leading to a disk- 
like accretion flow. For this reason, simulations to date 
have focused on disk-like accretion. However, it has been 
argued that for hot flows, in which the gas is near the 
galaxy virial temperature, the gas is supported by pres- 
sure against free infall and the flow may well be described 
by the spherical "Bondi" accretion model [||. In this 
so-called "cooling-flow model of quasar fueling", inter- 
galactic gas during the early stages of galaxy formation 
is accelerated toward the center of dark matter halos and 
is shock-heated to the virial temperature. This gas ac- 
cretes then onto the growing SMBH at nearly the Edding- 
ton rate. For a 10 6 M Q BH, the gas is expected to have a 
density n ~ 10 cm~ 3 and temperature T ~ 10 6 K — 10 s K 

Classical "Bondi accretion" refers to spherically sym- 
metric, steady-state accretion of an adiabatic gas onto 
a stationary star. The gas is assumed to be homoge- 
neous and at rest far from the star and flow adiabatically 
with a T-law equation of state (EOS). The problem was 
first studied in Newtonian gravitation for a point mass 
by Bondi [13], and later extended to accretion onto a 
BH in full general-relativity [H, [45^ ■ Accretion onto 
a star moving with constant velocity through a cloud 
which is asymptotically uniform and at rest was first 
studied qualitatively in Newtonian gravity by Hoylc and 
Lyttleton [5l| and later extended by Bondi and Hoyle 
[52|. The general-relativistic version for accretion onto 
a single BH has been studied via numerical simulations 
[53l [&4l [55I |56| . We refer to this scenario as the "Bondi- 
Hoyle-Lyttleton accretion" (BHL) problem. 

This paper is the first in a series of papers in which 
we will explore hydrodynamic gas flows around merging 
BHBH binaries. Here we focus on hot flows with zero net 
angular momentum and restrict our attention to simula- 
tions in which the binary is placed in a gas cloud which 
has constant density and temperature at infinity. We 
treat two cases: one in which the gas is asymptotically 
at rest (binary Bondi accretion), the other where the 
gas is moving with constant velocity with respect to the 
BHBH center of mass (binary BHL accretion). Bondi and 
BHL accretion onto single BHs represent classic prob- 
lems which are well understood. We seek to use this 
understanding as a foundation for tackling the problem 
of accretion onto a merging BHBH binary. The initial 
conditions of the gas during the merger of two supermas- 
sive BHs in realistic astrophysical environments is still 
an open question and subject to uncertainties in cosmo- 
logical structure and galaxy formation scenarios, and in 
our understanding of the formation history and role of 
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massive BHs in this process. The binary Bondi and BHL 
problems provide excellent settings in which to begin a 
rigorous probe of the binary accretion problem. 

The structure of the paper is as follows. In Sec. HT1 we 
discuss the unique computational challenge posed by the 
vast dynamical range characterizing our problem and we 
introduce our technique to tackle it. In Sees. IIIII and llVl 
we briefly outline the basic gravitational field and matter 
evolution equations and their specific implementation in 
our relativistic, hydrodynamic numerical scheme. Here 
we also provide an overview of our initial data, gauge 
conditions, and diagnostics. In Sec. IIV Dl we review 
code tests which were performed to validate our numer- 
ical scheme. In Sec. [V] we review the analytic Bondi 
solution for a single, stationary BH and demonstrate our 
ability to reproduce those results. We also compare with 
previous relativistic simulations of BHL accretion onto a 
single BH. In Sec. IVI1 we summarize the results of our 
binary BHBH merger simulations. In Sec. IVIII we sum- 
marize our findings and comment on future directions. 
We adopt geometrized units and set G = c = 1. 



II. COMPUTATIONAL CHALLENGE 

Simulating realistic accretion flows is extremely chal- 
lenging due to the enormous dynamic range in the rele- 
vant length scales defining the problem. One important 
length scale is set by the ADM mass M of the central 
gravitating source, R — M . This is the length scale at 
which relativistic effects become significant. Another im- 
portant length scale is the transonic radius, R s ~ M/a^, 
where aoo is the asymptotic sound speed in the gas cloud. 
This corresponds to the radius inside of which the accre- 
tion flow onto a stationary central mass M becomes su- 
personic. When the central gravitating object is moving 
supersonically with velocity Voo relative to the gas cloud, 
we follow [5l| and define the characteristic length scale, 
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Gas approaching the central object with impact param- 
eter less than ~ Rbhl will be captured. The two length 
scales R s and Rbhl can be combined into a single char- 
acteristic radius 14711. 
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within which gas is bound to M and will be accreted. 
Another important length scale for the BHBH problem is 
the binaray separation, d. When d <C R a , the accretion 
at radius r ^> d is basically BHL flow onto a central 
gravitating object of mass equal to the total ADM mass 
of the binary. When d 3> R a , the accretion in regions 
near each BH is again a BHL flow but for a source of 
mass equal to the mass of a single BH. 

The typical densities and temperatures for interstellar 
gas in a hot Bondi-like accretion flow ~ 10 cm 



and T ~ 10 6 K respectively [f|. Assuming that Voo ^5 a>oc, 
we find that R a ~ 10 6 M. With current computational 
resources, it is impossible to perform a relativistic sim- 
ulation that follows the complete binary inspiral from 
separation d ~ R a to d ~ M, assuming a realistic asymp- 
totic gas temperature. We address this issue by first per- 
forming "prototype" simulations in which we artificially 
increase the asymptotic temperatures in order to make 
the accretion radius R a much smaller. We then use these 
results and scaling to perform "realistic" simulations in 
which we treat realistic temperatures, but restrict our 
grid to domains much smaller than R a . The details of 
these approaches are given in Section [VTl 



III. BASIC EQUATIONS 

Throughout this paper, we use Latin indices to denote 
spatial components (1-3) and Greek indices to denote 
spacetime components (0-3). 

The basic gravitational field and relativistic hydrody- 
namics equations are discussed in [39L |57| where their 
numerical implementation is described and code tests are 
summarized. Here, we briefly sketch these equations and 
their implementation. 



A. Evolution of gravitational fields 

We write the spacetime metric in the standard 3 + 1 
form: 



ds 2 = 



2 dt 2 + 7 y (dx i + (3 l dt) (dx j + (3 j dt) . (3) 



where a, (3 l , and 7y are the lapse, shift, and spatial met- 
ric, respectively. The extrinsic curvature Kij is given by 



-2aKi 



(4) 



where Cp is the Lie derivative with respect to ff. We 
evolve 7,, and using the BSSN formulation [38l. [39t| - 
The fundamental variables for BSSN evolution are 



Hi 
K 



— ln[det(7y)] 

-id, 

e *7ij , 

ryV K ■ 



A 



\K l0 - - 7ij K) , 



(5) 

(6) 
(7) 

(8) 

(9) 



The evolution and constraint equations for these fields 
are summarized in [38l . [39(. We assume in this paper 
that the mass of the gas is negligible compared to the 
mass of the BHs, thus we do not include matter source 
terms in our metric evolution equations. 

Adding Kreiss-Oliger dissipation to the BSSN evolu- 
tion equations outside the BH can reduce high-frequency 
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numerical noise associated with AMR refinement inter- 
faces (43|. We use this technique here and have found 
it useful in reducing Hamiltonian and momentum con- 
straint violations. 

We adopt the standard puncture gauge conditions: 
an advective "1+log" slicing condition for lapse and a 
"Gamma- freezing" condition for shift [58] . Thus we have 

d a = -2aK, (10) 
do? = (3/4)5\ (11) 
d B l = d r - V B l , (12) 

where do = dt — ftdj. The r\ parameter is set to 0.5/M 
in all simulations. 

B. Evolution of hydrodynamic fields 

The fundamental matter variables are the rest-mass 
density po = nrriB, specific internal energy e, pressure 
P, and four- velocity u^. The stress-energy tensor for an 
ideal gas is given by 

T^u = pohu^Uv + Pg^u , 

where h = 1 + e + P/po is the specific enthalpy. We 
evolve the "conservative" variables p*, Si, and f. They 
are defined as 

p* = -^/ipon^ (13) 
Si = VlT^n^i (14) 
f = ^T^n^n" - p* . (15) 

Here n a = (a -1 , —a~ 1 (3 l ) is the timelike unit vector nor- 
mal to the t — constant time slices. Evolution equations 
are given by Eqs. (34), (36), and (38) of [13 ■ 

d t p* + dj(p*rf) = 0, (16) 

dtSi + djiay/^i) = ^a^T a Pd l9a0 ,(17) 

d t f + d l (a 2 ^T 0l -p*v l ) = s , (18) 

where 7 = det(jij) = e 12 ^ and v % = u l /u° is the fluid's 
3-velocity. The energy source term s is given by 

s = -a^T^V^ 

= a^[(T 00 (3 i p : > + 2T 0i ^' + T lj )K l3 

-(T oa f3 l +T 0l )d l a] . (19) 

C. Equation of State 

To complete the system of equations, we must specify 
an equation of state (EOS). While our code can handle 
any EOS of the form P = P(po, e), we adopt a standard 
r-law EOS, 

P = (r-l)p e. (20) 



We perform simulations with V — 4/3, 5/3, and 13/9. 
The choice of T = 4/3 is appropriate in the high- 
temperature limit in which both the electrons and 
baryons become relativistic (kT > m^c 2 ~ 10 13 if). The 
choice of r = 5/3 is appropriate in the opposite limit in 
which both electrons and baryons remain nonrelativis- 
tic (kT < m e c 2 - lO 10 ^). The choice of T = 13/9 is 
appropriate when the electrons are relativistic while the 
baryons remain nonrelativistic (m e c 2 ~ lO 10 !^ < kT < 
itibc 2 ~ 10 13 if)[il]. While it is not expected that tem- 
peratures in a realistic accretion flow will be high enough 
to make the matter behave as a T — 4/3 gas, we include 
it in our study as a limiting case (it may also be useful to 
model a radiation-dominated thermal g r = 4/3 

EOS). In general, T will also depend on the chemical 
composition of the gas. However, we show here that it 
is a good approximation to simply assume a pure gas of 
ionized hydrogen in assigning T for the matter-dominated 
regime of interest here. 

Take, for example, a fully ionized mixture of hydrogen 
and helium in which the electrons are relativistic while 
the nucleons remain nonrelativistic, and let X be the 
fractional abundance by number of hydrogen ions. The 
thermal energy density can be expressed as 

p e = X(3 + 3/2)nkT + (1 - X)(6 + 3/2)nkT . (21) 

Here, each relativistic electron contributes a factor of 
3nfcT, and each nonrelativistic nucleus contributes a fac- 
tor of (3/2)nfcT, where n = Po/tub is the baryon number 
density of the fluid. The pressure can be expressed as 

P = X2nkT + (1 - X)3nkT . (22) 

Combining Eq. ([201), Eq. ([IT]) and Eq. $Z2§ we find that 

_ (13/9)X + (7/3)(l-X) 

X+(5/3)(l-X) • [Z6) 

Adopting cosmological abundances [5^], we set X = 
0.92, which gives T = 1.439. Because this is close to 
T = 13/9 = 1.444, the value for X = 1, we henceforth 
set r = 13/9 for simplicity. In some cases, we allow 
for a temperature dependent transition from T = 5/3 
to r = 13/9 as the electrons become relativistic as they 
approach the BH horizons (see Sec. IV A 41 for details). 
Throughout this paper, we define temperature by 

P = 2nkT , (24) 

appropriate for pure ionized hydrogen. 

IV. NUMERICAL METHODS 

A. Initial Data 

For each simulation discussed in this paper, we con- 
sider the gas to be adiabatic (apart from shocks) with 
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uniform density and pressure at infinity. We take the gas 
either be asymptotically at rest (binary Bondi accretion) , 
or moving uniformly (binary BHL accretion) . We use the 
TWOPUNCTURES code [60| to construct the initial data for 
the binary BHBH metric. We choose the bare mass and 
momentum of the punctures according to [6l| in order to 
ensure that the BHBH binary orbit is initially close to 
quasicircular. 

We restrict our analysis to equal-mass, nonspinning 
BHs in this paper. However, we note that upon merger, 
the remnant settles down to a spinning black hole with 
ADM mass M f = 0.95AP Here M } denotes the final 
ADM mass, while M denotes the initial ADM mass of the 
binary. We also measure a final spin parameter a/Mf = 
0.69, in good agreement with the estimate of [62]. For 
this case, there is no recoil velocity. 

For all of our runs we use the analytic relativistic Bondi 
solution for accretion onto a stationary Schwarzschild BH 
as initial data. This solution is reviewed in Sec. [VJ To 
implement this data, we treat the binary as a single grav- 
itating object of mass M, where M is the ADM mass of 
the binary, and apply the analytic Bondi solution every- 
where outside a radius r$. We follow the method of [63[ 
and adjust the fluid parameters within the radius ro ac- 
cording to 



Po(r) = Po(r Q ) + 



dpp 
dr 



2r 



(25) 



This recipe ensures that the density and its first deriva- 
tive are continuous across ro. Eq. (|2"5|) is, of course, not 
the correct initial data for a quasistationary flow inside 
ro, but we find that the system settles into a quasista- 
tionary flow within several St ~ ro/a(r = ro), where a 
is the local sound speed. For all our BHBH runs, we set 
ro/M = 5.5, though we find that our results are not sen- 
sitive to this parameter. The fluid 4-velocity inside ro is 
set to be radially inward, with the magnitude set accord- 
ing to Eq. (f34|) . The fluid pressure is set according to a 
polytropic EOS, P — Kp^, with K = Poo/Poooj where 
Po,oo and Pqo are the rest-mass density and pressure at 
infinity, respectively. 

For our code tests in which we evolve only a single, 
stationary BH, our hydrodynamic initial data is con- 
structed in isotropic coordinates, which become singu- 
lar on the horizon at r = M/2. For these cases, we set 
ro/M — 0.6, and verify the finding of [63| that the fluid 
evolution quickly relaxes to the equilibrium Bondi solu- 
tion (see Sec. IV Bp . 



B. Evolution of metric and matter 

We evolve the BSSN field equations with fourth-order 
accurate, centered, finite-difference stencils, except on 
shift advection terms, where we use fourth-order accurate 
upwind stencils. We apply Sommerfeld outgoing wave 
boundary conditions to all BSSN fields. Our code is em- 
bedded in the Cactus parallelization framework [64| , and 



our fourth-order Runge-Kutta timestepping is managed 
by the MoL (Method of Lines) thorn, with a Courant- 
Friedrichs-Lewy (CFL) factor set to 0.5 in all BHBH 
simulations. We use the Carpet [65[ infrastructure to 
implement the moving-box adaptive mesh refinement. In 
all AMR simulations presented here, we use second-order 
temporal prolongation, coupled with fifth-order spatial 
prolongation. The apparent horizon (AH) of the BH is 
computed with the AHFinderDirect Cactus thorn [661 ] . 

We write the general relativistic hydrodynamics equa- 
tions in conservative form. They are evolved by a 
high-resolution shock-capturing (HRSC) technique [57j 
that employs the monotonized central (MC) reconstruc- 
tion scheme [67{ coupled to the Harten, Lax, and van 
Leer (HLL) approximate Riemann solver [6a ] . The 
adopted hydrodynamic scheme is second-order accurate 
for smooth flows, and first-order accurate when disconti- 
nuities (e.g. shocks) arise. Throughout the evolution, we 
impose limits on the pressure to prevent spurious heating 
and negative values of the internal energy e. Specifically, 
we require P min < P < P max , where 10 3 K p r 
and Pmi n = Kpq/2. Whenever P exceeds P max or drops 
below P min , we reset P to 

-^rnax 

or P m in, respectively. 
We check that this fix is applied only inside the apparent 
horizon, which is causally disconnected from the rest of 
the grid. 

At each timestep, we need to recover the "primitive 
variables" po, P, and v* from the "conservative" vari- 
ables p*, f, and Si. We perform the inversion as speci- 
fied in Eqs. (57)-(62) of [53], but with a slightly modified 
analytic quartic solver from the GNU Scientific Library 
that outputs only the real roots. We use the same tech- 
nique as in [69[ to ensure that the values of Si and f yield 

physically valid primitive variables, except we reset f to 
10 -io 

To, max (where To. max is the maximum value of f ini- 
tially) when either Si or f is unphysical [i.e., violate one 
of the inequalities (34) or (35) in [69|]]. The restrictions 
usually apply only to the region near the puncture inside 
the horizon. 



For all of our "prototype" calculations, we set our outer 
boundary at 328M and use 9 AMR refinement levels. 
The maximum resolution near each puncture is Sx/M = 
0.032. For our "realistic" calculations, we place our outer 
boundary at 164M and use 8 AMR refinement levels. For 
these cases, the highest resolution near each puncture is 
also Sx/M = 0.032. 

We model the emission of electromagnetic radiation by 
treating this radiation loss as a perturbation, and neglect 
its influence on the hydrodynamic flow as well as any 
deviation from adiabaticity. We also assume that the 
self-gravity of the gas can be ignored, and we do not 
include matter source terms in the evolution equations 
for the gravitational fields. 
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C. Diagnostics 

1. Rest-mass accretion rate 

An important diagnostic is the rest-mass flux through 
the horizon of each BH. In order to compute this quan- 
tity we begin by defining a function fit, x, y, z) such that 
fit, x, y, z) = on the 3d hypersurface which corresponds 
to the world tube of a BH horizon. Thus, we can write 

/ = V(x - x h (t)) 2 + (y- y h {t)) 2 + (z- z h {t)) 2 

-R(t,B,4>). (26) 

Here {xh(t), yh(t), Zh(t)} represents the coordinate posi- 
tion of the BH center at coordinate time t and R(t, 9, <fi) 
represents the coordinate distance from the BH center to 
horizon in the (9, (f>) direction. Here 9 and <fi are spherical 
polar coordinates with the origin set at the BH center. 
The rest-mass flux is then given by, 



M 



a v ^fp u fJ -d IJ ,fJd9d(l) , 



where J is the Jacobian 



J = 



d(t,f,6,4>) 


-l 




d(t,x,y,z) 




d(x,y,z) 



(27) 



(28) 



and where the integral is over each BH (apparent) hori- 
zon. 

The subtlety in this calculation arises from the fact 
that the position and shape of each horizon changes in 
time for a dynamical spacetime, so we cannot neglect the 
time derivatives of / in Eq. (|2T|) . For a full derivation of 
Eq. P?]). see Appendix 

We note that the M defined above is inherently gauge 
dependent. This is true even in a spacetime with a 
timelike Killing vector field. In our adopted gauge the 
sum of mass fluxes across the BH horizons is equal to 
the total mass flux through a large sphere measured 
by a distant observer in the cases in which the flow is 
(quasi)stationary. Quasi-stationary flow is realized dur- 
ing the BHBH inspiral phase, as well as after the binary 
merger once the remnant settles down. See Appendix [X] 
for a detail discussion on various issues involving M and 
gauge choices. 



In this paper, we instead adopt a crude method for esti- 
mating luminosities from our accreting gas, assuming it is 
an optically thin medium. Results should be interpreted 
as order of magnitude estimates only. Moreover, since we 
are primarily interested in enhancements and variations 
of the luminosity which may correlate with detectable 
gravitational wave signals, even crude estimates can pro- 
vide valuable information. 

We assume that the radiation propagates through an 
optically thin gas and we neglect the roles of radiation 
pressure and radiative cooling on the hydrodynamic evo- 
lution. Our solution is therefore valid only when the 
luminosity is below the Eddington limit and the radia- 
tive flux is below the thermal gas flux. To calibrate our 
crude approximation scheme we compare our results to 
the more exact calculation of optically thin emission from 
accretion onto a single, static Schwarzschild BH at rest 
in a gas cloud [4!|. We first review this solution below, 
then describe the approximations we make in this paper 
to simplify the calculation. 

The total luminosity received at infinity from gas 
accreting steadily onto a stationary Schwarzschild BH is 
given by , and accounts for Doppler and gravitational 
redshifts as well as photon capture from the BH. The 
result is 



dv 
dv 



dv 



(29) 



Here, dv^fdv accounts for Doppler and gravitational red- 
shifts, and is given by 



dv Q (1 - v 2 ) 1 / 2 



dv 1 — v cos 



7 (l-2Af/r) 1 / 2 



(30) 



where v is the proper velocity of the fluid as measured 
by a stationary observer, the specific luminosity L VQ is 
given by 



L Vn = 8tt 2 / r 2 dr 



2M 



cose* 

i. . 

I (1 — V COS 0') 



(1-v 2 ) 



and where 



|cos0* 



27/2M\ 2 /2Af \ 



1/2 



(31) 



(32) 



2. Luminosity 

In order to compute the observed electromagnetic ra- 
diation exactly, it would be necessary to employ a fully 
relativistic, radiative-transfer integrator suitable for a 
dynamical spacetime in 3 + 1 dimensions. While ad- 
vancements toward constructing such a scheme have been 
made in various a ppr oximations (e.g. [zO, [7l|, [72|, [ll, [zl, 
[ZE [ZE [Z3, H [ZIM HI, HI), this is an extremely chal- 
lenging task which has not yet been fully accomplished. 



Here the emissivity j„ is the specific emissivity (the en- 
ergy emitted isotropically per time per volume per fre- 
quency interval) in the comoving frame of the fluid. 

For dynamical systems containing inspiraling BHBH 
binaries, simple analytic expressions for the luminosity 
similar to Eq. (|3ip do not exist. In our simulations we 
employ a crude approximation for the luminosity. We 
simply compute 



L(t) 



j dV 



(33) 
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where dV is the proper 3-volume element in the fluid. 
We exclude from our integration points less than £r/j 
from each horizon. Here Th is the apparent horizon ra- 
dius, £ is a constant chosen so that for a single, isolated 
puncture, r — £rh corresponds to the surface of constant 
Schwarzschild radius f = 3M. We have made this choice 
because within this radius, 50% of the emitted radia- 
tionby a stationary source will be captured (see Eq. (24) 
in j49j|). For realistic temperatures (see Sec. IVI CI) , we 
find that the contribution to the luminosity is dominated 
by emission from near the binary. Our measured lumi- 
nosity is thus insensive to the outer limit of integration 
and we chose to integrate to the outer boundary of our 
integration domain for definitencss. For our high tem- 
perature prototype runs, however the gas has a nonneg- 
ligible contribution to the luminosity even far from the 
binary, thus the measurement is sensitive to the outer 
limit of the integration. We have chosen to integrate to 
r out — 44.7A/, which corresponds to r out « 2R a for our 
prototype BHBH Bondi runs. We have verified that for 
a single, stationary puncture, our luminosity estimate is 
within a factor of ~ 4 of the exact relativistic luminos- 
ity quoted in |49j , which we have separately verified with 
our code. Thus, our somewhat crude method of estimat- 
ing the luminosity is sufficient as an order of magnitude 
estimate. 

In our estimates of electromagnetic emission, we con- 
sider both thermal bremsstrahlung (free-free emission) as 
well as synchrotron emission. For bremsstrahlung emis- 
sion we consider both electron-electron and electron-ion 
processes. For synchrotron emission, we assume the pres- 
ence of a small-scale, turbulent B-field whose magnitude 
is approximated by setting P = (3Pm = (3B 2 /(8tt). We 
thus assume that the magnetic pressure is some fraction 
1/(3 of its equipartition value. Simulations of magne- 
tized accretion flows have demonstrated that the mag- 
netic fields do not typically reach their full equipartition 
value [83j . We have chosen /3 = 10 to account for this. We 
provide further details of the bremsstrahlung and syn- 
chrotron emissivities which we use in Appendix [B] 



D. Code Tests 



All of the above tests and simulations were performed 
on grids with uniform spacing. In some of the simula- 
tions, we utilized the multiple-transition fisheye transfor- 
mation I87| so that a uniform computational grid spacing 



corresponds to physical coordinates with spatially vary- 
ing resolution. Recently, we have modified our code so 
that we can use the moving-box AMR infrastructure pro- 
vided by Carpet [65| . To test our new code, we have per- 
formed shock-tube tests and 3+1 simulations of linear 
gravitational waves, single stationary and boosted punc- 
ture BHs, puncture BHBH binaries, and rapidly and dif- 
ferentially rotating relativistic stars. Our AMR code has 
also been used to perform simulations of BHNS mergers 

All of our 3+1 AMR code tests were performed as- 
suming equatorial symmetry (i.e., symmetry about the 
z = orbital plane), which we assume in all evolutions 
presented in this paper. We have also checked that our 
AMR code is able to accurately reproduce the analytic 
Bondi solution, and we have shown good agreement with 
the results of [53| for BHL accretion. Results of the latter 
two tests are summarized in Sec. IV Bl 



V. ACCRETION ONTO A SINGLE BH 

A. Analytic Relativistic Bondi Solution 

Steady state, adiabatic, spherically symmetric accre- 
tion onto point masses was first considered by Bondi 
for Newtonian gravitation (47[. This work was later ex- 
tended to handle accretion onto BHs in general relativ- 
ity [H, Hi| [H, Hl| . A thorough discussion of the relativis- 
tic solution may be found in Appendix G of [5(| • Here we 
briefly summarize this solution. We consider spherically 
symmetric, steady-state accretion onto a Schwarzschild 
BH of mass M. We assume that the BH is at rest in an 
infinite gas cloud which has rest-mass density p^, pres- 
sure Poo, and fluid 4- velocity u 1 — at infinity. We as- 
sume that the gas is adiabatic with an adiabatic index T. 
We can solve the equations of relativistic hydrodynamics 
to derive an exact solution for the accretion flow. 



Our general relativistic hydrodynamic code has 
been thoroughly tested by passing a robust suite of 
tests. These tests include maintaining stable rotating 
stars in stationary equilibrium, reproducing the exact 
Oppenhcimer-Snyder solution for collapse to a BH, and 
reproducing analytic solutions for relativistic shocks and 
spherical Bondi accretion onto isolated BHs [57j . Our 
code has also been used to simulate the collapse of 
very massive, rotating stars to black hole s 18411 ; merging 
BHBH binaries H, BHNS binaries [H and rela- 
tivistic hydrodynamic matter in the presence of puncture 
black holes [63j . Recently, our code has been generalized 
to incorporate (optically thick) radiation transport and 
its feedback on fluids in dynamical spacetimes [82|. 



1. Key Equations 

For simplicity, we will derive the solution in 
Schwarzschild coordinates, then transform to isotropic 
coordinates. We begin by recasting the steady-state rel- 
ativistic continuity and Euler equations into an integral 
form, 

4irp a uf 2 = M = const. (34) 
h2 ( X ~ ~~T + = h lo= const. (35) 

Here, f and u denote the radius and radial fluid 4- velocity 
in Schwarzschild coordinates. 
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For adiabatic flow, the T-law EOS implies the poly- 
tropic relation P = Kp\, with K =const. Thus 



p = Kpi = (r - i) Po 



h-l 

r 



For an adiabatic gas, the speed of sound is, 



h 1 / 2 



dP\ 1/2 (TP 
dpo) \poh 



1/2 



(36) 



(37) 



This leads to the following relations between the en- 
thalpy, the speed of sound, and the temperature 



a 2 = (r-i 



h-l 



kT 



P m B r - 1 



2n 



(h-l). 



(38) 
(39) 



Any solution satisfying Eqs. (|34|) and (|35| that maintains 
the causality constraint a 2 < 1 must contain a transonic 
radius outside the event horizon At the transonic 

radius the radial velocity ii s and the speed of sound are 
given by 



M 
2R S 



M/2R S 



1 - 3M/2R S 



(40) 
(41) 



The transonic radius sets the length scale at which the 
fluid parameters begin to deviate from their asymptotic 
values during inward flow. We derive the transonic radius 
by plugging Eq. (|4i"j) into Eq. (G.30) from Appendix G 
of [50[. This leads to a polynomial equation which we 
solve using a numerical root solver. In the Newtonian 
limit, in which <C 1 and M/R s <C 1, we recover the 
expressions, 



R„ = 



5 - 3r M 

— S£ lsr<5/3 

3 M 

4 o m 



5/3 . 



(42) 



Knowing the transonic radius, the accretion rate is 
given by 



M = 4np s u s R 2 s = 4ir\M 2 p oa a 



-3 



(43) 



where A = A(r, q,oo) is a parameter of order unity, which 
in the Newtonian limit is given by 



A ='5 



j \ (r+i)/2(r-i) / 5 _ 3r 



-(5-3r)/2(r-i) 



■ (44) 



In the relativistic domain, we must solve a cubic equation 
for to find A (see [5(| for details). Given values for p .oo 
and Too, we can determine M and h^ using Eqs. (f38|) . 
(1551) . and (|25]). We can then use Eqs. ^ and ([5S). along 
with the EOS, Eq. (f3"6")l . to obtain the complete analytic 
solution. 



2. Transformation to isotropic coordinates 

In order to provide initial data for our simulations, we 
transform the solution to isotropic coordinates, the coor- 
dinate system adopted for our metric (puncture) initial 
data. The relation between Schwarzschild radius f and 
isotropic radius r outside the horizon is given by 



M 
l + Yr 



f — M + y/f(f - 2M ) 



Since the flow is purely radial, we have 



(45) 
(46) 



f dr 
df 
u 



(1 - M/2r)(l + M/2r) ' 



(47) 



where u r is the radial component of the 4-velocity. The 
time component of the 4-velocity, u , remains unchanged 
by the transformation. We find the 3-velocity v — v r — 
u r /u from the normalization condition u a u a = —1, or 



a 2 {u») 2 + e 4 *{u°) 2 v 2 



-1 



(48) 



The above transformation is needed only when assign- 
ing fluid initial data for r > rQ. Inside tq we set the 
initial data according to Eq. (f25|) . We choose r = 0.6M, 
which is outside the horizon, to avoid the initial coordi- 
nate singularity of u° at the horizon. There is no singu- 
larity during the evolution as the adopted 1+log slicing 
condition is horizon penetrating. 



3. Moving Bondi solutions 

We also wish to consider the problem of adiabatic, 
steady-state accretion onto a BH which moves with con- 
stant velocity Voo through a gas which is uniform at 
infinity. This problem was first studied in Newtonian 
gravitation in [5l| and [52| . and is often referred to as 
"Bondi-Hoyle-Lyttleton (BHL) accretion" . The problem 
has also been studied numerically in Newtonian theor y by 
[5J, |55j , and for BHs in full general relativity by 
For cases in which the gas is not at rest asymptotically, 
we construct initial data following the method described 
in 63], in which we "boost" the stationary Bondi solu- 
tion by taking the initial density at any coordinate point 
to be approximately the stationary isotropic Bondi solu- 
tion, and computing the initial velocity field employing a 
Lorentz transformation for the velocities. Denoting the 
stationary Bondi flow velocity by v l , the moving Bondi 
flow by v l , and the "boost" velocity by Voo, we set 



ct 



1 + 



(49) 
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FIG. 1: Semi-analytic profiles of density, temperature and 
fluid 4-velocity for spherical, adiabatic Bondi accretion onto 
a Schwarzschild BH. The solid line gives po/po,oo, the dot- 
ted line shows T/Toa, and the dashed line shows the fluid 
4-velocity in Schwarzschild coordinates, u. The adiabatic in- 
dex is given by T = F* according to Eq. (|53[) . The asymptotic 



4- Effective adiabatic index 

In some cases, we consider a gas cloud which the elec- 
trons have nonrelativistic temperatures at infinity, but 
achieve relativistic temperatures near the binary. In this 
case we follow [4!| and define an "effective adiabatic in- 
dex" r* according to 



r* 



f 5/3 
I 13/9 



kT/m e c 2 < 2/3 
kT/m e c 2 > 2/3 



(53) 



Thereby replacing the actual continuous transition by a 
simpler discrete transition. This transition occurs at a 
Schwarzschild radius [501 



r. 
M 



9 m p 
40 m e 



tt « ~ 400 



(54) 



temperature is Too = 10 K. The horizon is labeled Ru- 



We solve for using Eqs. ([gg]). (j3"5)) . (|3"9"|) . and the fact 
that fcT(r*) = (2/3)m e c 2 . Using the continuity of po 
and P at r = r*, we obtain the full equilibrium flow 
solution. In practice, r* is outside the outer boundary 
of the computational grid in our simulations, so we still 
implement a constant T. However, taking into account 
this transition alters our initial data significantly, since 
the outer T = 5/3 region drives a < u up to r = r*, 
increasing the gas temperature near the black hole when 
compared to gas in which T = 13/9 everywhere. Analytic 
profiles of density, temperature, and fluid velocity for this 
equation of state are plotted in Fig. Q] 



B. Relativistic Bondi test 



Cl 



76(1 + Kx>f) 



(50) 
(51) 



where 7;, = (1 — V^,) -1 / 2 and c; = a/ip 2 . Here a is the 
lapse, and ip = e^, with defined in Eq. ([5]). For any 
hydrodynamic quantity g(x,y, z), we compute 



g{x',y',z') = g(j b x,y,z) , 



(52) 



which accounts for the coordinate transformation be- 
tween the x a and x a frame on the t' = timeslice. The 
above method of computing initial data is only strictly 
valid asymptotically. This is because we use the Bowen- 
York initial data for the moving BH metric, which are not 
obtained by boosting the stationary BH metric. How- 
ever, we find that deviations propagate away quickly, 
leaving a stationary flow. In Sec. IV Bl we describe a sim- 
ulation of Bondi accretion which we perform in a moving 
reference frame. In this case, we use the same technique 
described above to construct hydrodynamic initial data. 



We test our code's ability to accurately simulate hy- 
drodynamic accretion onto a moving puncture. This is 
particularly important in an AMR code in which matter 
crosses moving refinement zone boundaries. To test our 
code, we simulated spherical Bondi accretion in a frame 
in which both the BH and gas cloud are moving with the 
same velocity. We consider a BH moving with velocity 
Vbh = P x /M = 0.37, initially situated at the origin. 
Here P x is the momentum of the BH as measured by a 
stationary coordinate observer at infinity. We consider a 
gas with adiabatic index T = 13/9 with asymptotic sound 
speed ale = 0.022 also moving with velocity v = Vbh at 
infinity. With this choice the transonic radius is given 
by areal radius R a = 45. 5M in the comoving frame. We 
evolve for a duration t = 713. 6M, which is approximately 
equal to two free-fall times at the accretion radius. By 
this time, the BH has moved to a coordinate location of 
x = 247 M. 

To assess the agreement with the analytic solution, we 
compare invariant quantities [63j |. Two such invariants 
are the fluid rest-mass density, po, and the the rate of 
change of the fluid rest-mass density, as measured by a 
comoving observer, dpo/dr. We choose a set of 6 areal 
radii in the comoving frame. At each radius, we com- 
pute po and dpo/dr analytically. By construction, con- 
tours of po and dpo/dr are spherical and coincide in 




FIG. 2: Snapshots showing spherical Bondi accretion onto a 
single BH simulated in a frame in which the BH and gas move 
with velocity Vbh = 0.37. Contours of constant density po 
(solid black lines) and constant dpo/dr (dotted red lines) are 
shown. Their overlap indicates agreement with the analytic 
Bondi solution. The adiabatic index is set to V = 13/9 and 
a~o = 0.148. The green dashed lines show the boundaries of 
the innermost AMR refinement levels. Arrows denote velocity 
vectors. 



this frame. However, because po and dpo/dr are both 
coordinate-independent quantities, their contours must 
continue to coincide in any coordinate frame. We use this 
fact to check that our numerical simulation, performed 
in a frame in which the BH and gas are moving, matches 
the analytic solution (see Fig. [2J . 



C. Relativistic Bondi-Hoyle-Lyttleton test 

We will simulate cases of accretion onto binaries in 
which the BHBH center of mass moves relative to the 
asymptotic gas cloud. Here we verify that we can accu- 
rately simulate "BHL" accretion of a single puncture BH 
moving at constant velocity through a gas cloud. While 
there is no analytic solution for this problem, we com- 
pare our findings with results of previous work. We con- 
sider a test case in which a single BH puncture is placed 
in a cloud with asymptotic sound speed a<x, = 0.1 and 
r = 4/3, with the puncture moving supersonically with 
speed Voo = 0.25 relative to the gas cloud. We perform 
our simulation in a frame in which the puncture is at rest 
and the asymptotic fluid velocity is set to V^. Figure [3] 
shows a snapshot of the stationary-state density and ve- 
locity profile from our simulation. We measure M and 



FIG. 3: Snapshots of density contours for the BHL ac- 
cretion code test. The adiabatic index is T — 4/3, the 
sound speed is <2oo = 0.1. Density contours are chosen at 
Po = po.oolO ' 25 -' (j = 1, 2, ....12). Contours of highest density 
are very near the BH. Arrows represent velocity vectors. 

compare with the results of [53||. In order to facilitate 
this comparison, we define a canonical unit of rest-mass 
accretion flux to be 

Mean- {vl+alaf/2 , (55) 

We find M / M cll n = 2.6, which is slightly smaller than 
the value 3.0 reported in [53j. The most likely sources of 
the small discrepancy are the outer boundary condition 
(we use the stationary spherical Bondi solution with a 
Lorentz boost, whereas 15311 use constant asymptotic val- 
ues), and the fact that [53| impose an approximate inner 
boundary condition outside the horizon. Other differ- 
ences between these two codes are that we use a 3+1 code 
with AMR, whereas [53| employ an axisymmetric code. 
Also, our outer boundary is placed at r max /M = 820 
for this test, whereas [53j place the outer boundary at 
r max /M = 140. 

VI. ACCRETION ONTO A BHBH BINARY 

As discussed in Sec. HH the transonic radius for Bondi 
accretion with a realistic asymptotic temperature of 
10 6 K is R a ~ 10 6 M. It is beyond the capability of cur- 
rent 3+1 GR simulations to evolve the binary from initial 
separation d > R a all the way to merger, while resolving 
a BH horizon of size ~ M . The range of length scales is 
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TABLE I: Parameters for BHBH simulations 



run run type case Voa/a x V TqJJQ a x x RajM n d/Af 



PAl 
PA2 
PA3 


prototype 


Bondi 


0.0 


0.0 


13/9 
4/3 
5/3 


8.73 x 10 iu 
9.62 x 10 10 
7.44 x 10 10 


0.148 


22.7 


40, 20, 14, 10 


PBl 
PB2 
PB3 




BHL 


0.1 


0.7 


13/9 
4/3 
5/3 


8.73 x 10 iu 
9.62 x 10 10 
7.44 x 10 10 


0.148 


15.6 


40,10 


PCI 
PC2 
PC3 




BHL 


0.4 


2.7 


13/9 
4/3 
5/3 


8.73 x 10 1U 
9.62 x 10 10 
7.44 x 10 10 


0.148 


2.74 


40,10 


RAl 
RA2 


realistic 


Bondi 


0.0 


0.0 


r* 

5/3 


1 x 10 b 


5.53 x 10~ 4 


3.27 x 10 b 
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R a — (M/2)/(a^o + V^) 3 ' 2 is the accretion radius for a single BH of mass M/2. 
'd — initial binary separation; simulations beginning at 10M proceed to merger. 



TABLE II: Parameters for BHBH simulations 



run flow characteristics^ emission characteristics' 







Tjnax 1 Too 






f f t max j r 


l max 

nv is 


tt j max 1 j 
^syn 1 -^35 


U v max 






PAl 


631 


18.3 


1.2 


3.7 














PA2 


1743 


14.3 


1.5 


3.8 














PA3 


164 


30.1 


1.6 


2.6 














PBl 


593 


18.9 


1.2 


4.0 














PB2 


1638 


15.2 


1.7 


3.3 






Not Physically Relevant 








PB3 


153 


30.8 


1.9 


2.5 














PCI 


99.2 


28.3 


6.0 


0.7 














PC2 


158 


32.2 


14.0 


0.7 














PC3 


46.6 


37.7 


6.5 


0.8 














RAl 
RA2 


1.4 x 10 10 
3.63 x 10 9 


1.8 x 10 6 
2.7 x 10 6 


4.7 
13.0 


3.9 
3.3 


300 
400 


150 MeV 
230 MeV 


3 x lO 8 /^ 1 80/(1 + z) 
IxlO 8 ^ 1 100/(1 + 2) 


n\ /2 T b 


3/4^-1/2 
-3/4^-1/2 


eV 
eV 



^M a c 2 = 4.6 x 10 40 A 5 / 3 Hi T 5 ~ 3/2 Mi erg s" 1 is the accretion rate onto a single BH of mass M/2 
undergoing stationary, spherical Bondi accretion. 
t1 "L 35 = 10 35 n{ T 6 - 3 M| erg s" 1 

§ "max" label refers to the characteristic value at the moment of maximum luminosity. 

§§ "max" label refers to the maximum values of the luminosities during the merger and the 

characteristic frequencies at these times. 



too large and the total coalescence time too long for such 
a task. We approach this issue by performing two types 
of simulations. We perform simulations of binaries merg- 
ing in "realistic" gas clouds with asymptotic temperature 
Too = 10 6 K, following only the last phase of the merger 
in which the binary separation satisfies d <C R a - Our 
focus here is on identifying observable electromagnetic 
signals generated by the time-dependent shock heating 
caused by the binary motion. We also perform "proto- 
type" calculations with artificially high temperatures and 
sound speeds in order to study how the accretion flow 
changes as the binary transitions between the following 
regimes during the inspiral: 

1. "widely separated regime", in which d > R a 

2. "moderately separated regime", in which d ~ R a 

3. "closely separated regime" , in which d < R a 



We set the asymptotic temperature ~ 10 11 if in the 
"prototype" calculations to decrease the accretion radius 
R a to a value closer to d so that we can explore all three 
regimes by combining the results of several numerical 
simulations. All of our BHBH simulations are summa- 
rized in Table |U Important results are summarized in 
Table [II] The initial BHs in the binary are all equal- 
mass and nonspinning. 



A. Scaling 

For a given asymptotic gas temperature or sound 
speed aoc, our solution for binary Bondi flow can be 
scaled to arbitrary density = po.oo/ws (neglecting 
self-gravity of the gas) and black hole mass M. Hence 
a single simulation with an arbitrary rirx, and M suf- 
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fices to determine the solution for any other or M. 
Thus, for example, the accretion rate is proportional to 
««, and M 2 (see, e.g. Eq. (|43| ). while the 4- velocity u a 
as a function of coordinates x a /M are independent of 
rioo and M . If the asymptotic sound speed is sufficiently 
low that doo <C 1, then the solution can also be scaled 
to arbitrary or T^. However, once approaches 
the speed of light and the transonic radius approaches 
the horizon, scaling with or breaks down. This 
behavior is already evident from the relativistic Bondi 
solution for accretion onto a single BH. 

The emergent electromagnetic luminosity and radia- 
tion spectrum also exhibit simple scaling. The form 
of the scaling relations depend on the temperature and 
density dependence of the adopted emissivities and will 
be described later in Sec. IVI CI when we treat realistic 
asymptotic temperatures and sound speeds. 



B. Prototype Cases 

1. Binary Bondi Accretion 

The "prototype" calculations provide valuable qualita- 
tive insight into the evolution of the accretion flow as the 
binary separation decreases, passing through the three 
regimes identified above. For the case of "binary Bondi 
accretion" , in which the gas at infinity is at rest relative 
to the binary center of mass, we obtain a sequence of 
"snapshots" for different binary separations. Each snap- 
shot is generated by evolving the binary long enough to 
allow the gas to relax to a quasistationary flow, but not 
long enough for the separation d to change significantly. 
By sandwiching together snapshots we can trace the evo- 
lution in accretion rate M and luminosity L as the bi- 
nary transitions from regime 1 to regime 2. The final 
transition from regime 2 to regime 3 is captured by a 
single simulation which follows the binary through inspi- 
ral and merger. We perform these simulations for cases 
PA1 (r = 13/9), PA2 (r = 4/3), and PA3 (r = 5/3) 
in order to study the effect of the EOS on the evolution 
(see Table HJ . Numerical results for all cases are given in 
Table [TTJ The time variations in the accretion rate and 
total luminosity for these cases are shown in Figs. 2H6] 

The evolution in M can be understood as follows. De- 
fine M a to be the accretion rate onto a single, isolated 
BH of mass M/2. Then it follows that the total accre- 
tion rate onto two infinitely separated BHs of mass M/2 
will be given by M — 2M a . However, late in the inspi- 
ral, when the separation satisfies d <C R a , the binary can 
be treated as a single gravitating object. From Eq. (|43|) 
we see the accretion rate is proportional to M 2 , so we 
expect that the accretion rate will approach M = 4M a . 
During the final stage of the inspiral, mass-energy is ra- 
diated away in the form of gravitational waves. Thus we 
expect that the final post-merger accretion rate will be 
approximately given by M = 4M a (l — SM/M) 2 . In our 
simulations, we find 8M/M ss 0.05, consistent with the 



d / M 

40 20 10 8 




t / M 

d / M 

40 20 10 8 




FIG. 4: Time evolution of M and SL for binary Bondi ac- 
cretion for the prototype case with F = 13/9. Time t/M is 
measured relative to the time at which the merger occurs. 
M a and SL a are the accretion rate and luminosity enhance- 
ment over the background for a single isolated black hole with 
mass equal to the initial ADM mass of the binary. The left- 
hand box shows values from "snapshots" in regimes 1 and 
2. The right-hand box shows results for the final inspiral 
and merger from an initial separation of d = 10 M (regime 3). 
Solid lines denote numerical data, dashed lines denote extrap- 
olated data. Solid dots correspond to profile plots highlighted 
in Fig. [7] while open circles show expected values at t = ±oo. 
The asymptotic sound speed is set at doo = 0.148, for which 
Ra/M = 22.7. 



values reported in [62j. We observe the expected behav- 
ior for both our PA1 (see Fig. @} and PA2 (see Fig. [5]) 
runs. For PA3, the accretion rate never reaches AM a be- 
fore the merger. We attribute this to the fact that for 
r = 5/3, gas is more efficiently heated, allowing P ~ pov 2 
and a ~ v, even in the absence of shocks. When shocks 
do form, the flow is much more easily disrupted as the 
kinetic energy of the flow does not dominate the thermal 
energy, contrary to cases PA1 and PA2. Thus, some mat- 
ter is swept away from the vicinity of the binary, causing 
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FIG. 5: Same as Fig.H but for T = 4/3. 



FIG. 6: Same as Fig.H but for F = 5/3. 



a lowering of the accretion rate (see Fig. [6|). We note, 
however, that after the merger when the shocks have a 
chance to dissipate, the accretion rate does settle to its 
expected value, taking mass loss into account. 

We also plot the luminosity enhancement due to 
bremsstrahlung and synchrotron emission in Figs. 
Because the high-temperature homogeneous background 
gas in our prototype simulations has an intrinsic, non- 
negligible emissivity, we subtract it from the total lumi- 
nosity measured. We define SL = L — Lb g , where Lb g 
is the background luminosity which would be present in 
our computational domain for a homogeneous gas cloud 
of density po.oc and temperature X^, with no BH present. 
We normalize SL by SL a , which we define as the lumi- 
nosity above background which would be present for a 
single, isolated BH of mass M/2. If we take the limit 
in which the binary separation d — ► oo, we expect that 
SL/SL a — > 2. This limiting value is indicated by the 
open circle at t/M = — oo in Figs.[4]-[6] We also calculate 
the expected value of 5L/SL a for a single BH of mass 
M — SM and plot it for reference with an open circle at 
t/M — +oo. We see that for each r the luminosity en- 



hancement increases by several orders of magnitude over 
the course of the inspiral. While the numerical value of 
this variation is not physically meaningful due to the un- 
realistic temperatures used in these "prototype" calcula- 
tions, this behavior provides strong qualitative evidence 
of a significant enhancement in luminosity that can be 
expected to accompany such an inspiral. Mergers in re- 
alistic clouds, yielding realistic luminosities whose values 
are physically meaningful, will be treated in Sec. IVI CI 

Figures [7] and [5] show snapshots of density and tem- 
perature contours for cases PAl and PA3. We do not 
show snapshots for the PA2 case because they look very 
similar to the PAl case. We can see that in the early 
phases of the inspiral, the accretion flow resembles two 
independent spherical Bondi flows. As the separation de- 
creases and becomes comparable to the transonic radius, 
the orbital velocity of each BH becomes comparable to 
the sound speed. Within R a shocks begin to form, which 
grow in strength until the merger. It is the heating from 
these shocks which contributes to the dramatic increase 
in the luminosity observed. Note that the final accre- 
tion flow near the BH is not spherically symmetric due 
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to the spin of the merged BH. The BH spin does not 
significantly change the final accretion rate predicted by 
Eq. (|43|) . as this quantity is determined by gas parame- 
ters at r ~ R a ^S> M , where the effect of the BH spin is 
negligible. This result conforms with the findings of [89j . 

In order to highlight the role that shock heating plays 
during phases 2 and 3 of the merger, we also present 
contours of K/K^ for our PA1 case (see Fig. [9|). Here 
K = P/pq and Koo is the value of K at infinity. Be- 
cause K/Koo = 1 everywhere for adiabatic flow in the 
absence of shocks, this quantity serves as a useful tracer 
for the amount of shock heating which is taking place. 
The quantity K — K(s), where s is the gas entropy, and 
is constant in the absence of shocks; shock heating yields 
K/Koo > 1 (see Appendix B of [86]). As expected, we 
see that KjK x increases steeply near the shock front. 
For each snapshot, we compute the maximum value of 
K/Koo outside the horizon. We find that K max j K^ ini- 
tially increases as the separation decreases and shocks 
become stronger, as expected. This trend terminates in 
the very late stage of the merger, when d <~ 5M ~ disco 
[9fj| . This is likely due to the fact that kinetic energy dis- 
sipated as heat is confined to a small region and is being 
quickly consumed by the BHs at this stage. After the 
merger, the gas relaxes to laminar spherical Bondi flow 
and K/Koo returns to unity everywhere. 



2. Binary Bondi-Hoyle-Lyttleton Accretion 

In order to investigate additional electromagnetic sig- 
natures which may be present due to the motion of the 
binary relative to the cloud, we have performed a series 
of "prototype" simulations of BHL accretion onto merg- 
ing binaries. We consider both subsonic cases (Voo/aoo = 
0.7) and supersonic cases (Voo/o-oo — 2.7). We again con- 
sider T = 13/9 (PB1 and PCI), V = 4/3 (PB2 and PC2), 
and r = 5/3 (PB3 and PC3) in order to assess the in- 
fluence of the EOS on the flow. We perform both wide 
separation runs (d — 40M) to produce snapshots that 
are quasistationary in the corotating frame of the binary, 
as well as close separation runs (d = 10M) in which we 
evolve the binary from inspiral to merger. 

For our (asymptotically) subsonic cases, the departures 
from the binary Bondi case are subtle. Figures ITTJ1 and [TT1 
show snapshots of density and velocity field for cases PB1 
and PB3. Snapshots for case PB2 are similar to PB1 and 
hence are now shown here. At wide separation d — AQM , 
we see an asymmetry in the accretion flow for cases PB1 
and PB2, as shocks develop around one BH as it moves 
against the flow of the gas, but not around the other as 
it moves in the same direction as the flow. This phe- 
nomenon continues up to the merger. At this separation 
the orbital Keplerian velocity is Vk ~ 0.08. Thus, these 
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FIG. 10: Same as Fig. [3 but for subsonic BHL accretion with F — 13/9 and Voo = 0.1. The asymptotic velocity Voo is in the 
+x direction. 
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FIG. 12: Time evolution of M and 8L for binary BHL inspi- 
rals with T = 13/9, = 0.148, and Vao = 0.1. The initial 
binary separation is d = 10M and the BHs evolve to merger. 
M a and SL a are the accretion rate and luminosity enhance- 
ment over the background for a single isolated black hole with 
mass equal to the initial ADM mass of the binary. Dashed 
lines in M plot represent accretion rates onto individual BHs, 
the solid line represents the total accretion rate. Dots rep- 
resent the times highlighted in the last three snapshots in 
Fig. [TO] 



shocks are formed when one BH moves supersonically 
against the flow of the gas and (V/ c + V 00 )/a > 1. The rea- 
son that this behavior is not seen in case PB3 is that the 
gas is adiabatically heated more efficiently for T = 5/3 
as it flows toward the BHs, causing the sound speed of 
the gas to be greater near each BH. We find that for 
r = 5/3, (Vfc + Vao) J a < 1 near the black hole, prevent- 
ing any shocks from forming. The post-merger accretion 
flow exhibits some departure from spherical symmetry 
due to BH spin and ^ 0, but all shocks dissipate. In 
Figs. [12UT41 we once again observe an increase in accre- 
tion rate and luminosity over the course of the merger. 
We have also plotted in these figures the individual ac- 
cretion rates onto each BH in order to demonstrate the 
effect of the binary motion relative to the wind. 

For our supersonic cases (PC1-PC3), the departure 
from the binary Bondi case is more dramatic. When 
the BHs are widely separated (d > R a ), a bow shock 



FIG. 13: Same as Fig.ll2l but for F = 4/3. 

forms around each individual BH, as seen in Figs. [TBI and 
[TBI Late in the inspiral, when d < R a , these shocks 
merge and form a single bow shock surrounding the bi- 
nary, which persists after the merg er. The final flow re- 
sembles the solution found by [531 ] for steady accretion 
onto a moving BH with Voo/onx = 2.5, although in our 
case the remnant is spinning. During the inspiral, the 
modulation in M and SL is more pronounced than in the 
subsonic case (see Figs. [T7HT9|) . 



C. Realistic Binary Bondi accretion 

While the high temperature prototype runs give us 
valuable insight into the general nature and different 
phases of the accretion flows onto inspiralling BHBH bi- 
nary systems in gas clouds, we require simulations with 
more realistic gas temperatures in order to identify ob- 
servational signatures. Accordingly, we have performed 
simulations of the binary Bondi problem for a gas cloud 
with asymptotic density tIqq = 10 cm -3 and temperature 
T = 10 6 K. This choice is consistent with the proposed 
"cooling flow model of quasar fueling" describing the hot 
interstellar gas found in galaxies [a EH- As mentioned 
previously, computational limitations demand that our 
outer boundary be placed inside the transonic radius for 
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these simulations. As a result, we focus only on the final 
phase of the inspiral and merger, in which d <C R a . We 
use both a T = T* (case R Al) an d r = 5/3 (case RA2) 
EOS. As explained in Sec. IV A 41 we still set T = 13/9 
in the computational domain for the r = T* case since 
the r = 5/3 and transition region is outside our com- 
putational domain. However, the initial hydrodynamic 
profile of the flow is very different from a pure T — 13/9 
EOS with the same asymptotic temperature. We focus 
here on the BHBH Bondi problem only (V^ = 0), and 
postpone a study of binary BHL accretion in a realistic 
gas cloud for a later analysis. As in the prototype calcu- 
lations, we find evidence for a strong enhancement in the 
luminosity due to shock heating of the gas (see Fig. [20] 
and Fig. l2~Tj) . For these runs, we plot the luminosities in 
cgs units. 

We note that for case RAl, strong shocks form, but are 
confined to the immediate vicinity of the binary. This is 
because the inward gas flow onto the binary is strongly 
supersonic, making it difficult for shocks to propagate 
outward (see Fig. [22]) . For case RA2, on the other hand, 
the radial component of the 4- velocity u ~ a everywhere, 
making it much easier for shocks to spread outward, as 
seen in Fig. |2"3"1 



D. Scaling and Detectability 

All of our quoted results for the accretion rate M 
are normalized by the value for a single black hole of 
mass M/2 undergoing stationary, spherical Bondi accre- 
tion. Using Eq. ([43]) we see that this quantity scales with 
asymptotic temperature and density according to 



M„ 



■? = 4.6 x 10 40 A 5/3 rii T 6 3/2 M? 



erg s 



(56) 



A 



5/3 



A(r,a oc )/A(5/3,a 

~|6] 



fix 

\6 1 



where we define 

noo/lOcm" 3 , T 6 = T oo /10 6 K, and M 6 = M/1O 6 M . 
Recall that in the Newtonian limit, when R a ^> M, A 
depends only on T. 

It is also possible to derive simple scaling relations for 
the luminosities. In each of our simulations, most of the 
electromagnetic luminosity is generated near the horizons 
of the BHs, as the temperature and density both rise 
sharply when approaching the horizon, and achieve their 
maximum values there. By examining Eqs. (]B1[) - (]B5|) 
(ignoring the logarithmic terms), and Eq. (|B18[) . we see 
that in the high temperature limit (kT > m e c 2 ), which 
applies near the horizon in all of our simulations, the 
bremsstrahlung and synchrotron emissivities depend on 
the density and temperature according to, 



Iff « n\T h 



Qsyn n h 



(57) 
(58) 



Here and Th refer to the density and temperature at 
the horizon, respectively. Integrating Eqs. ([57]) and (|58[) . 
we find 



L fJ « J dVq ff cx n 2 h T h M 3 



(59) 



L S yn ps / dVq syn cx nj^/^Af 3 , (60) 



where we have assumed that the radiation is generated 
near the horizon and hence the radiative volume scales 
as M 3 . To estimate nh, we use the fact that r/, R a , 
and so the fluid 4-velocity is approximately given by its 
free-fall value 



1/2 



Substituting this into Eq. 



Here we have used a 2 



3/2 



and Eq. 
( 

\ itibc 2 



(61) 



we find that 



-3/2 



(62) 



« TP/p a = 2TkT/m B c 2 . This 
scaling should remain reasonably accurate, even in the 
presence of strong shocks, as the shocks cause density 
enhancements of < (r + l)/(r - 1) (see Sec. 89 of (HI), 
which is of order unity. 

The temperature at the horizon T; l7 on the other hand, 
will be strongly affected by the presence of shocks. The 
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FIG. 15: Same as Fig. [7] but for subsonic BHL accretion with F — 13/9 and Voo = 0.4. The asymptotic velocity Voo is in the 
+x direction. 
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FIG. 16: Same as Fig. [7] but for subsonic BHL accretion with T — 5/3 and Voo = 0.4. The asymptotic velocity Voo is in the 
+x direction. 
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0.4. 



FIG. 18: Same as Fig. [T5J but for V = 4/3, = 0.4. 



gas in any shocked region will be heated to kT ~ m B v 2 . 
Since v ^ c near the horizon, shock heating guaran- 
tees that kTh ^ rase 2 ~ 10 13 K independant of the gas 
temperature at infinity T^. This result has important 
consequences for the scaling of the maximum luminosi- 
ties. Once the shock-heated gas is accreted following the 
merger, Th drops below this value for all T < 5/3 (see 
e.g. Fig Q]) and the luminosity settles down to a value 
below the maximum. 

Even in the absence of shocks, kTh ^ mgc 2 for 
r = 5/3. This is because for this particular EOS, an 
appreciable fraction of the gravitational potential energy 
is converted into thermal energy (both scale as r^ 1 inside 
R a : kT ~ Mm B /r) For r = T*, it can be shown 

that the temperature at the horizon for spherical Bondi 
flow is approximately given by [50] 



kT h 



2/3 



m B 



1/3 



m B c 2 = 0.013 m B c 2 (63) 



Thus, we find that the temperature at the horizon, Th is 
independent of gas parameters at infinity for both T = 
5/3 and T = T* . However, for T = const < 5/3, T h does 
exhibit scaling with Too in the absence of shocks. In this 



case, P — Kp , so we find 

r-i 



\rriBC 2 



-3(r-i)/2 



(64) 



We can now apply these results to see how the lu- 
minosity scales in different phases of the inspiral. Dur- 
ing both the very early pre-merger phase of the inspiral, 
when d > R a , and during the post-merger phase after the 
fluid has settled to a quasiequilibrium state, there are no 
shocks present. We can therefore use the above results 
to see that in these phases, 



n\ T 6 - 3 M| , r = 5/3 or V = V* , 

Lff i n\ T^ T+1)/2 Mt , T = const < 5/3 , (65) 



oc 



-9(r-i)/ 2/ ,-i,. 3 



T = 5/3 or T = F 



n\ Tf J{L - 1>/z p^ l M'i , T - const < 5/3 



(66) 



Here (3\ = (3/10. During the late premerger phase 
when d ~ M <C R a shocks will be present and Th will no 
longer depend on Too for any EOS. In this case, we find 

L ff <xnl T 6 " 3 M 6 3 , L syn cx n\ T^^M 3 . (67) 

Using the above scaling relations along with the results 
of our realistic temperature simulations, we find that the 
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FIG. 19: Same as Fig. [H but for F = 5/3, V x 



0.4. 



peak luminosity shortly before the merger of an equal 
mass BHBH binary in a gas cloud for case RAl (r = T*) 
is given by 

3 x 10 37 n\ T 6 - 3 M| erg s" 1 , (68) 
3 x 10 43 n\ T^P^Ml erg s^ 1 . (69) 



j max 

L ff 

t max 
syn 



Similarly, the peak luminosity for case RA2 (F = 5/3) is 
given by 

4 x 10 37 n\ T 6 3 Ml erg s" 1 , (70) 
4 x 10 43 n\ Tq 3 /3^MI erg s^ 1 . (71) 



j max 

syn 



At the late post-merger phase, the fluid relaxes to a 
stationary flow. The scaling relations j65|) and (|66| hold. 
Combining these scaling relations and our simulation re- 
sults, we find that during the post-merger phase 



L ff w 3 x 10 35 n\ T 6 - 3 M| erg s" 1 , 
L syn « 8 x 10 38 n\ T f r 3 /3r 1 Af 3 erg s" 

for case RAl, and 

Lff k. 3 x 10 36 n\ T 6 _3 M| erg s" 1 , 
L syn « 2 x 10 41 n\ T^P^M* erg s" 



(72) 
(73) 

(74) 
(75) 



for case RA2. We note that in each case, the total lumi- 
nosity is dominated by the synchrotron emission. 



FIG. 20: Plots showing time evolution of M and L. Here 
time is measured relative to the time at which the merger 
occurs. Asymptotic temperature is T = l(fK. Adiabatic 



index given by V ■ 
M 6 = M/10 6 M Q . 



,/lOcm-^, T 6 = Too/10 6 K, 



In each of our calculations, we have ignored the effects 
of radiative cooling on the gas dynamics. We can esti- 
mate the error induced by this by comparing the rate of 
thermal energy transport, Eth, to the luminosity. Here 
we define 



Eth = Me = M 



po(r - 1) 



(76) 



Since the luminosity is dominated by emission near 
the horizon, we are primarily concerned about the region 
near the horizon. Using Eqs. 03]), §2$, G2) and 

(fT3")) . we find that for case RAl, at late times after the 
merger when the flow has reached equilibrium, 



L 



sy n 



L 



E 



^ ~ 0.1 ni T 6 3/2 p^Me 



(77) 



Thus, we see that in these regimes, it is a good approxi- 
mation to neglect the effects of radiative cooling for our 
canonical model parameters. At the moment of maxi- 
mum luminosity, shortly before merger, we find that 



40 m Tg" 372 ^" 1 M, 



E th 



(78) 
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FIG. 21: Same as Fig.ECTl but for F = 5/3. 



Thus during the final stages of the merger, the validity 
of our assumption of adiabatic flow begins to break down 
for canonical parameters. In future work, we will address 
this by including cooling terms in our gas evolution to 
account for energy losses due to radiation. 

We also note that in order for us to be able to neglect 
radiation pressure in the momentum equation, we require 
that the luminosity be small compared to the Eddington 
luminosity. We find that for case RA1, 



t max I j max 



L 



Edd 



0.2 n\ T 6 - 3 /3f l M { 



(79) 



which suggests that radiation pressure may begin to play 
a role for parameters close to our canonical choices. 

In calculating the luminosity, we have assumed that 
the gas is optically thin. We can verify this assumption 
by estimating the optical depth. For the gas parameters 
chosen for this study ( Tirxj — 10 cm 3 and Tqq — 10^K), 
we find that the dominant source of opacity is electron 
scattering. We estimate that the optical depth for elec- 
tron scattering of synchrotron photons is 



n h a T R 

-,-3/2 



10~ 3 niT fi 3/2 Af 6 



(80) 



where ~ 10 niT 6 is the density at the horizon, 
(Tt = 0.67 x 10~ 24 cm 2 is the Thomson scattering cross- 



section, and R ~ 10 11 Mg cm is the characteristic size of 
the emission region. Thus, our assumption of an optically 
thin gas is valid. 

To estimate the characteristic frequencies at which 
the emission occurs we again note that the maximum 
emission comes from near the horizons and compute the 
characteristic frequency produced in this region. For 
bremsstrahlung emission, the characteristic observed fre- 
quency of the emission is given by hv ~ kTh/(l + z) 
for a source at redshift z. We measure the maximum 
temperature near the horizon for our case RAl, and find 
that at the moment of maximum luminosity in the late 
pre-merger phase, 



l. ..max 

hv Sf 



150 MeV 
l + z 



(RAl) 



(81) 



After the merger, in the quasistationary phase, we find 
that the characteristic frequency drops to 



hv ff 



10 MeV 
l + z 



(RAl) 



(82) 



Following the same procedure for case RA2, we find that 
at the moment of maximum luminosity in the late pre- 
merger phase, 



nU ff 



230 MeV 



(RA2) 



l + z 

and in the post-merger phase, the frequency drops to, 
50 MeV 



(83) 



1 



(RA2) 



(84) 



Thus, we see that the bremsstrahlung emission will 
be predominantly in 7-rays, in agreement with [49| . 
Given the bremsstrahlung luminosity calculated above, 
we estimate that the flux from this emission will be 
~ 10~ 21 erg cm~ 2 s _1 for a source at z = 1. Unfortu- 
nately, it is unlikely that this emission is strong enough 
to be detectable. 

We can estimate the characteristic frequency of the 
synchrotron emission by noting that Eq. (|B10|) is maxi- 
mized when xm ~ 1.09. For case RAl, this corresponds 
to an observed frequency 



U v max 
"'"syn 



1.09 3ehB ( kT 



l + z 47TTO e c V rn e c 2 
80 



(85) 



1 



n l/2 r _3/4g-l/a e y (RA1) (8g) 



during the late pre-merger phase at the moment of maxi- 
mum luminosity. In the post-merger phase, the frequency 
drops to, 



hp 



syn 



0.008 
l + z 



For case RA2, we find the characteristic synchrotron fre- 
quency to be 



h v max 



100 1/2 ^-3/4^-1/2 



T 6 Pi eV ( RA2 ) 
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during the late pre-merger phase at the moment of max- 
imum luminosity, and 

hvsyn = n\' 2 T~ 3/4 p; 1/2 eV (RA2) (89) 

during the post-merger phase. 

This corresponds to infrared and visible radiation. For 
a binary at z = 1 with the luminosity calculated above, 
this source has an apparent magnitude of m = 24 and 
should be observable by the proposed LSST instrument 
[92| . Our simulations follow the late stage of the inspiral 
in which the binary separation decreases from d — 10M 
to merger. For a 10 6 M Q binary, this corresponds to a 
timescale of At ~ 1.3 hrs during which the radiation 
should achieve peak values. 

We note that all the scalings derived above will not ap- 
ply when Voo 3> 

&00 1 but can be derived in a similar fash- 
ion. For a realistic gas with T > 1Q 6 K (a^ w lOOkm/s), 
this regime will might never be realized, so we neglect 
here. 



VII. SUMMARY AND DISCUSSION 

In this paper we have performed a set of fully gen- 
eral relativistic simulations of BHBH binary mergers in 
gaseous environments, using our relativistic hydrody- 
namics code with AMR capability. Our focus has been 
on demonstrating a mechanism for an observable elec- 
tromagnetic signal which may accompany the gravita- 
tional wave signal from a black hole merger. We have re- 
stricted our attention to gas clouds which are asymptoti- 
cally uniform and either at rest (BHBH Bondi accretion) , 
or moving (BHBH BHL accretion) with respect to the 
binary center of mass. We have performed "prototype" 
high-temperature (T^ ~ 10 11 K) simulations in order 
to gain a qualitative understanding of the different flow 
regimes characterizing such systems. In these simula- 
tions, the accretion radius R a is placed within the compu- 
tational domain, which allows us to sample all the various 
regimes, studying how the luminosity and accretion rate 
change when the binary separation d goes from d > R a 
to d < R a . We have also performed "realistic" simula- 
tions with astrophysically plausible asymptotic temper- 
atures (T!^ = 10 6 K) to identify observable electromag- 
netic signals. For each simulation, we have calculated 
the time-varying rest-mass accretion rate, as well as the 
electromagnetic luminosity due to bremsstrahlung and 
synchrotron emission. We also derive scaling relations 
for the luminosity for the realistic cases, enabling our re- 
sults to be used for different asymptotic gas parameters 
and BH masses. 

In each case, we find evidence for a time-varying elec- 
tromagnetic signature accompanying the BHBH binary 
merger. For our "realistic" temperature simulations, 
we find that during the final inspiral from separation 
d = 10M to merger, there is an enhancement in the 
total luminosity of ~ 1 — 2 orders of magnitude. This 



is followed by a sharp drop in the luminosity by a fac- 
tor of ~ 2 — 5 orders of magnitude immediately following 
the merger. In each case, we find that the luminosity is 
dominated by the synchrotron component, and that this 
signature should be detectable by the proposed LSST in- 
strument. 

We note that the nature of the electromagnetic emis- 
sion may be altered significantly when one considers dif- 
ferent accretion geometries, such as disk accretion, for ex- 
ample. We intend to study such systems in future work. 

By solving the BHBH binary analogs of the classic 
Bondi and BHL accretion problems we have sought to lay 
a natural and rigorous foundation for our future simula- 
tions involving BHBH mergers in gaseous environments. 
Our focus here was to perform "prototype" simulations to 
identify the different regimes characterizing the gas flow 
in these classic scenarios and then to perform a "real- 
istic" simulation for a merger occuring in an astrophys- 
ically plausible gaseous setting. For the latter case we 
also determined the scaling behavior of the gas dynam- 
ical parameters, emitted luminosities and characteristic 
frequencies so that they could be used to apply our re- 
sults to other environments. We noted that for asymp- 
totic gas parameters not too different from the ones that 
characterize the interstellar medium in recent models of 
merging galaxies (rico ~ 10 cm -3 and ~ 10 6 K) some 
of the assumptions that go into our calculation may be 
breaking down, like the assumptions of adiabaticity, sub- 
Eddington luminosities and an optically thin medium. 
These are complications we intend to address in future 
work. We also hope to compare our findings to other re- 
cent investigations performed simultaneously with ours 
(e.g. [13), although these treat a different scenario and 
adopt very different initial data, making a direct com- 
parison difficult. Finally, we hope to use the simulations 
reported here as point of comparison with future simula- 
tions involving BHBH mergers in ambient gaseous disks. 
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APPENDIX A: DERIVATION OF M 
EXPRESSION 

Consider a 3D hypersurface F given by f(t,x,y,z) = 
constant in a spacetime diagram (see Fig. l24[) . The inter- 
section of F with a t = constant time slice is a 2-surface 
S(t). Below we will identify S(t) to be the apparent hori- 
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S(t +At) 



t +At„ 




FIG. 24: Spacetime diagram depicting the hypersurfaces rel- 
evant to calculating the mass accretion rate. The 2D hyper- 
surface S(to) and S(to + At) (white regions) represent BH 
horizons on neighboring time slices. They are enclosed by 
spacelike 3D hypersurfaces Et and £t +At (shaded region) 
on these slices. 



zon of a black hole at time t. Let L be a 3D hypersurface 
which is a worldtube enclosing F. Let us further define 
S t to be the 3D region on the time slice t between the 
surface S(t) and S7(i) , the 2D cross section of L on the 
time slice t. We can imagine a 4-volume V to be the re- 
gion bounded by the hypersurfaces S to , T,t +st, L, and F. 
Consider a 4-vcctor flux j 7 ' with a vanishing divergence, 
V^j 71 = 0. At a given time t, define the function q(t) 
according to 



?(*) = 



J 



^ 3 n„ 



(Al) 



From Gauss's law, 

o = / v^dv = [ rd 3 n, 

JV JdV 
= q(t + St) - q{t ) 



^d 3 n, 



(A2) 



To evaluate the integral on F, it is convenient to 
introduce a coordinate system (t,f,a,b), where a — 
a(t, x, y, z) and b = b(t, x, y, z) are two other coordinates. 
We may write 



^ 3 n„ 







3! J 




3! j 


f j f £ftab 


-J 





dx v dx p dx° 



= - J V^fdnfJ dtdadb , (A3) 

where g' is the determinant of the metric in the (t, /, a, b) 
coordinate system, g is the determinant in the (t, x, y, z) 
coordinate system, J is the Jacobian 



J 



d(tj,a,b) 


-l 


d(f,a,b) 


-l 


d(t,x,y,z) 




d(x,y, z) 





(A4) 



and jf — j^d^f follows from the usual transformation 
formula for a vector field between the (t, /, a, b) and 
(t, x, y, z) coordinate systems. Similarly, if L is given by 
l(t, x, y, z) — 0, then 



^d 3 U^ = - 



-gj^d^Ui dtdadb 



(A5) 



where Ji = \d(l, a, b)/d(x, y, z)\ 1 . Taking the limit 8t 
0, we obtain 



dq 
lit 



-T F + J~ L 



where 




gj^dfilJi dadb . 



(A6) 

(A7) 
(A8) 



Consider a fluid accreting onto a BH. Let F be the 
BH horizon world tube. At any given time t, consider 
the fluid in the region E t between the BH (apparent) 
horizon S(t) and a distant 2-surface ft(t) surrounding the 
BH. Let L be the 3D hypersurface formed by stacking f2 
with time. The continuity equation gives V ^(pou^) = 0. 
Setting j M = pou^, we have 



,d 3 x 



M (t) (A9) 



is the rest mass bounded by the surface f2 and the hori- 



zon, where p* 



-gpQU . Hence we have 



dMp 
dt 



(A10) 



This equation states that the rate of change of the rest 
mass is equal to the amount of rest mass flowing into O 
per unit time (JFl) minus the amount of rest mass flowing 
into the horizon per unit time (Tf)- Hence we define the 
rest-mass accretion rate onto the BH according to 



M = T F 



a^/jpouV-dftfJdddcf) 



(AH) 



which is the expression given in Eq. (|27|) . Here we have 
used the identity y/—g — oty/j, and choose a and b to be 
the spherical angular coordinates 8 and cf> with the origin 
at the BH center. 

It is apparent from the definition that in general M 
depends on how the spacetime is sliced near the horizon. 
However, in some cases M may be time independent. 
Consider the cases where the fluid's mass is negligible 
compared to the BH's mass and there exists a timelike 
Killing vector £ = d/d\ in the vicinity of a BH. This 
Killing vector could be the time Killing vector describing 
a stationary BH, or a helical Killing vector which approx- 
imates the BHBH spacetime in the inspiral phase. One 
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might choose to measure M using a coordinate system in 
which t = A (at least locally). In this case, 



M = M\ = — / y/^QxPoU 3 
IF 



(A12) 



where g\ is the determinant of the spacetime metric 
in the (A, /, 0, (j>) coordinate system. Suppose the flow 
of the fluid also achieves a stationary state near the 
horizon in which d\(pou^) = everywhere on F, i.e. 
pou? — pQuf(6,4>) on F. Since dxgx = 0, M\ is a con- 
stant independent of the coordinate time A. On the other 
hand, in a coordinate system in which A is not the time 
coordinate, we have 



M = - ^/^g 7 P ou f d6dct > = - / (ey'V^Po^dOd^, 

J F J F 

(A13) 

where g' = gx/i^) 2 is the determinant of the spacetime 
metric in the (t, f,9,<f) coordinate system, and is the 
time component of the Killing vector £. Note that both 
g\ and pou? are still time independent on F, but M is 
time dependent if £* is time dependent. The accretion 
rate M is time independent only if a gauge is chosen 
so that d t £* = on F. Furthermore, if £* is constant 
everywhere on F, then M = Mx/^, where is the 
value of £* on F. As an example, consider a stationary 
accretion flow onto a Kerr BH. In th e boost ed Kerr-Schild 
coordinates, we have £* = 7& 



everywhere in 



the spacetime, where Vf, is the boost velocity. Hence we 
have M = Mx/lb in the boosted Kerr-Schild coordinates. 

We point out that in a numerical simulation, even if 
a Killing vector exists, the adopted gauge (i.e. time slic- 
ing) may not correspond to the gauge in which d/dt is 
the Killing vector. However, in some situations there ex- 
ists a gauge in which a (quasi)stationary flow is expected. 
Such situations include the Bondi accretion onto a BHBH 
binary in the inspiring phase and the Bondi accretion 
onto a single BH following the binary merger. In these 
situations, the mass accretion rate onto a distant, fixed 
surface SI (i.e. Fl) is time independent for any gauge 
choices that give rise to a spacetime that is asymptoti- 
cally Minkowsky. A "well-behaved" gauge should give a 
Fp (or the sum of two Fp^s in the BHBH case) equal to 
Fl (when averaged over time); otherwise, Eq. (|A10|) im- 
plies that there will be an accumulation (if Fp < Fl) or 
depletion (if Fp > Fl) of rest mass in the interior of SI as 
a result of a pure gauge effect. In our numerical simula- 
tions, we do not see such a gauge effect. We compute Fl 
on spherical surfaces of various radii. We find that after 
the flow reaches a (quasi)stationary state, Fl is slowly 
changing with time in the binary inspirling phase and is 
approximately time independent after merger. The com- 
puted fluxes at various radii are also the same. Moreover, 
the sum of the computed fluxes at the BH horizons agree 
with the value Fl, indicating that our adopted puncture 
gauge conditions are well-behaved. 



APPENDIX B: EMISSIVITIES 

a. Bremsstrahlung emissivity 

In order to estimate the electromagnetic emission due 
to bremsstrahlung, we use the following expressions for 
electron-ion, and electron-electron cooling rates given in 

M 



Iff = + q ee 
,8tt 



(Bl) 



q el = n 2 — {a f r 2 e c) {m e c 2 )F a {9) ergs cm s (B2) 
q ee = n 2 (a f rlc)(m e c 2 )F ee (9) ergs cm~ 3 s _1 (B3) 

Here, a/ = e 2 /hc is the fine structure constant, r e = 
e 2 /m e c 2 is the classical electron radius, 9 = kT/m e c 2 , 
n = po/ tub is the baryon number density, and 
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1/2 



4( — l (1 + 1.78 1 1 - 34 ) 9 <1 

/-,■(*, = [ y) 

— [ln(l. 1230 + 0.48) + 1.5] 6 > 1 
2tt 



(B4) 



F ee {9) 



2 i_ (44 _ 3^3/2 



97T 1 / 



(1 + 1.10 + 6 2 - 1.25 5/2 ) 0<l( B5 ) 



240(ln 1.1230+1.28) 



b. Synchrotron emissivity 



> 1 



We use the estimates for synchrotron cooling rates 
given by [94| : 



where 



Anne 2 v T (x M \ _ 3 _ 1 _ a 

— = 1 ergs cm s Hz , 

V3cK 2 (l/9) Vsin0/ 

(B6) 

= cyclotron frequency (B7) 

(B8) 



vq = 

X M = 



eB 

2v 
Zv 9 2 



From [951 ] . we get the following approximation for I(xm)- 



I(x m ) = 2.561 1 + — 



1.92 0.9977 



/3 1 2/3 
\ X M X M 

and the angle averaged version, 



exp(-1.8899x] v { 3 ) , 
(B9) 



„ . 4.0505 / 0.40 0.5316 \ , 1/3n 
I'(x M ) = 1 + — + exp(-1.8899x;f ) 

X M V X M X M ) 

(B10) 

We integrate over frequency to get the total cooling 
rate. 
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Let us denote, 



Here we have used 



4.0505 



Anne 



VscK 2 (i/{e)) 



01 



a 2 = 1.8899 



(BH) 

(B12) 
(B13) 



Thus, we see, 

> 

q v $dv 



A 



n Jo 



5/6 



0.4 0.5316 



■ M 



1/4 

5 M 



x exp 



(-a 2 :r] v ,{ 3 ) dx M 



SAC 



Where 



C 



r(ll/2) 0.4T(19/4) 0.5316 T(4) 



11/2 



19/4 



Thus, wc find 



1/2 I 

C M I 

(B14) 

2.151889 
(B15) 



' 108r? 2 e 4 r 5 

<Zs = I q v , S dv = 4.0505 C ^ ^ , ns n 9 . (B16) 



P = $P M = /?f- . 

sir 



(B17) 



As discussed in Sec. IIV1 we have chosen = 10 for our 
simulations based on simulations of magnetized accretion 
flows which have demonstrated that the magnetic fields 
do not typically reach their full equipartition value [83| . 
We further assume that the hydrodynamic flow is tur- 
bulent due to the magneto-rotational instability (MRI), 
causing the frozen-in field lines to be tangled, and ran- 
domly oriented [88]. This justifies our use of the angle 
averaged function in Eq. (|B10[) . Note that for a; <C 1, 
K 2 (x) w 2 jx 1 , so that for sufficiently large temperatures 
we may approximate: 



q Vt8 dv = 4.0505 C 



54n 2 e 4 c 9 3 
^3 (3m e c 2 



(B18) 
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